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Abstract 

The cosmic ray interaction event generator sibyll is widely used in extensive air shower sim- 
ulations. We describe in detail the properties of sibyll 2.1 and the differences with the original 
version 1.7. The major structural improvements are the possibility to have multiple soft inter- 
actions, introduction of new parton density functions, and an improved treatment of diffraction. 
sibyll 2.1 gives better agreement with hxed target and collider data, especially for the inelastic 
cross sections and multiplicities of secondary particles. Shortcomings and suggestions for future 
improvements are also discussed. 
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I. INTRODUCTION 



Cosmic ray interactions in the atmosphere can be regarded as high energy fixed target 
collisions involving heavy particles. Because of their low intensity, cosmic rays with energies 
above 10 15 eV can only be studied indirectly through the extensive air showers (EAS) they 
initiate in the atmosphere. The analysis of EAS relies on air shower Monte Carlo simulations 
which uses hadronic interaction models. At higher energies, where the cosmic ray energy 
is beyond the reach of man-made accelerators, hadronic interaction properties have to be 
extrapolated. The difficulties in the extrapolation are augmented by the fact that, while 
the forward region contains most of the energetics and is important for shower development, 
most of the accelerator measurements are made in the central region. 

The event generator sibyll [1] is intended for air shower cascade simulations. It is a 
relatively simple model that is able to reproduce many features of hadronic interactions in 
fixed target and collider experiments, sibyll is based on the dual parton model (DPM) [2-4] , 
the Lund Monte Carlo algorithms [5, 6], and the minijet model [7-10]. The hard interaction 
cross section is calculated according to the minijet model. For hadron-nucleus interactions, 
the interaction probability for each nucleon inside the nucleus is calculated based on the 
impact parameter distribution. The total interaction cross section is calculated using the 
Glauber scattering theory [11]. For a nucleus-nucleus interaction the semisuperposition 
model [12] is used to determine the point of first interaction for the nucleons of the projectile 
nucleus. The fragmentation region is emphasized as appropriate for air shower simulations. 
Versions 1.6 and 1.7 of sibyll have been released and used since the early 1990s. The only 
difference between the two is that version 1.7 can have neutral pion interactions, which is 
important only for air showers above 10 19 eV because at lower energy all neutral pions decay 
before they interact. 

Several shortcomings of sibyll 1.6 and 1.7 have been identified over the years, such as 
(1) the total proton-proton cross section calculated with the parton structure functions rose 
faster than what the experimental measurements indicate; (2) multiplicity fluctuations and 
average charged particle multiplicity are too small at high energy; (3) diffractive events did 
not agree well enough with the available data sets. For these reasons the event generator 
was modified and has been available as sibyll 2.1 [13] since 1999. 

The most important changes in version 2.1 are in the description of soft interactions 
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and diffraction dissociation. In order to allow multiple soft interactions, the eikonal for the 
soft interaction is described using Regge theory, whereas in version 1.7 the eikonal for the 
soft interactions was energy independent and had the same b dependence (b is the impact 
parameter) as used for hard interactions. While in version 1.7 the cross section for diffraction 
dissociation is parametrized independently of the eikonal model, a two-channel eikonal model 
based on the Good- Walker model [14, 15] is used in sibyll 2.1. In addition, low- and high- 
mass diffraction dissociation are treated separately in the new version. However, it should be 
kept in mind that diffraction dissociation is still not satisfactorily understood. The parton 
structure functions have been updated to agree with the HERA results. Other parameters 
were retuned with updated values as well. The multiple soft interaction and new parton 
densities give larger multiplicity at high energies and better agreement with data. The 
multiplicity distribution has been improved by implementing better the effect of diffraction 
dissociation. 

The aim of this paper is to describe the current 2.1 version of sibyll to make a reference 
of the implemented physics models and ideas available. We will outline the overall structure 
and improvements made, within details of the soft interactions and diffraction dissociation. 
We compare sibyll with fixed target and collider data, and we show how it performs in air 
shower simulations. Finally, we list some remaining shortcomings of sibyll 2.1 and outline 
how they can be improved. 



II. HADRON-HADRON INTERACTION 



A. Basic DPM picture 

sibyll 2.1 retains the DPM picture. In the DPM picture, a nucleon consists of a quark (q, 
color triplet) and diquark (qq, color antitriplet). Soft gluons are exchanged in an interaction 
and the color field gets reorganized. The projectile quark (diquark) combines with the target 
diquark (quark) to form two strings. Each string fragments separately following the Lund 
string fragmentation model [6]. 

The fractional energy x of the quark f q (x) is chosen from a distribution of 

where a = 3.0 and \i = 0.35 GeV is the effective quark mass. The diquark energy fraction 
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is then f qq {x) = 1 — f q {x). If particles 1, 2 collide to form strings a and b, the energy and 
momentum of the strings are as follows 
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To fragment the string, a g-g pair or qq-qq pair is generated at one of the randomly chosen 
ends of the string. The new flavor combines with the existing one to form a hadron, and 
the remaining (anti)fiavor becomes the new end. A primordial p? of equal magnitude and 
opposite signs is assigned to the pairs, with a Gaussian distribution where the mean is energy 
dependent 
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with p = 0.3 (u, d), 0.45 (s), 0.6 (qq) and ^Js being the cm. energy of the hadron-hadron 
interaction. The energy fraction of each new particle follows the Lund fragmentation function 
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where a = 0.5 and b = 0.8, rriT = \J m? + p\ is the transverse mass, and z is the fraction 
of the new particle energy with respect to the parent quark or diquark. The fragmentation 
process continues until the remaining string mass is less than a "threshold mass." The 
threshold mass here is defined as the quark masses of the string ends plus the quark/diquark 
pair mass plus (1.1 ± 0.2) GeV. The string finishes the fragmentation by forming two final 
hadrons. 



B. Hard interactions and minijets 

Already in the range of collider energies, at \fs ~ 100 GeV, the original DPM picture of 
just two strings cannot explain what is observed, namely: 

1. high multiplicity; 

2. increase of mean px] 

3. high p T jets; 

4. rise of central rapidity density. 
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These new features can be interpreted as the emergence of hard interactions which become 
prominent as energy increases, in the form of minijets. Minijets are described with pertur- 
bative QCD, where partons from the colliding hadrons experience hard scattering. Minijets 
have a transverse momentum larger than some momentum transfer scale p™ in ^> Aqcd, 
where perturbative calculation holds, but smaller than a typically reconstructed collider jet. 
The minijet formalism described below is based on Refs. [1, 16] with modifications made in 
the new version. 

Minijets are perceived as part of the hard interaction described by perturbative QCD. 
The cross section is calculated within the QCD-improved parton model in leading order is 



VQCD(s,PT m ) = K j ' dxi J dx 2 J dp T 
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where f a ,i{ x i-,Q 2 ) an d fb,j(x2,Q 2 ) are the parton distribution functions of parton i (j) in 
particle a (b). The transverse momentum of the scattered partons is denoted by p?- The 
calculation is done for four light flavors. Higher order corrections are accounted for by setting 
the factor K=2 and the factorization scale Q 2 = p\. sibyll 1.7 used parton densities of 
Ref. [17] (EHLQ), where the gluon density is extrapolated as g(x) ~ 1/x at small x. Data 
from HERA [18, 19] suggest a steeper increase at low x. sibyll 2.1 uses parton densities 
of Ref. [20, 21] (GRV) which scales the gluon density as l/x 1+A with A = 0.3 - 0.4. As in 
the previous version, Eq. (6) has been calculated separately and is included in the code in 
tabular form. 

The change in the low-a; region affects the minijet cross section substantially at high en- 
ergies. The cross section cannot rise without limit at high energies [22, 23]. If the number of 
gluons times the transverse resolution scale of hard interaction (~ 1/p^) becomes compara- 
ble to the proton size, nonlinear effects, possibly saturation, cannot be neglected. Another 
factor to take into account is the use of collinear factorization approximation in calculating 
minijet cross sections, where the transverse momenta of the incoming partons (i, j) should 
always be smaller than the transverse momenta of the scattered partons (k, I). This approx- 
imation is used to sum the parton densities and only the leading term ln(p^) is considered. 
The collinear factorization approximation breaks down for ln(l/x) 3> ln(p|,). The ln(l/x) 
term becomes important at high energies and needs to be taken into account [24]. In order 



to restrict the calculation of the minijet cross section to the region of phase space where the 
QCD-improved parton model is expected to be reliable, the following transverse momentum 
cut is applied 

p^(s) = p° T + Aexp c^/ln(s/GeV 2 ) , (7) 

where p^ = 1 GeV, A = 0.065 GeV, c = 0.9. This parametrization follows from the geometric 
saturation condition [23] 

^-.xg{x,p%) < nR 2 p , (8) 
Pt 

where a s is the strong coupling constant, g(x,p^) the gluon density, and R p is the ef- 
fective radius of a proton in transverse space. The scale Q 2 = p\ is assumed. In the 
limit ln(l/x), ln(Q 2 /A 2 ) — > oo (double leading-logarithmic approximation) the steeply ris- 
ing gluon density g(x, Q 2 ) can be written 



xg(x, Q ) ~ exp 
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with A being the QCD renormalization scale and rif is the number of quark flavors. The 
functional form of Eq. (7) follows from inserting Eq. (9) in Eq. (8), however, the parameters 
in Eq. (7) cannot be derived directly from first principles. 

SIBYLL 1.7 had a constant cutoff at p™ m = a/5 GeV. Using modern parton 
density parametrizations one cannot obtain a satisfactory description of the proton- 
proton/antiproton cross sections with a transverse momentum cutoff that is energy- 
independent. The py m values for both versions are shown in Fig. 1. 

The minijet cross section quickly rises to exceed the total cross section. This is interpreted 
as the collision forming more than one minijet. The average number of hard interactions 
n hard occurring at energy s and at impact parameter b is [9] 

^hardO, s) = A(b)a QC D(s) , (10) 

where A(b) is the profile function for the hadron-hadron collision. The baryon and meson 
profile functions follow those of Refs. [10, 25] and are given in Appendix B. Following the 
convention given in Ref. [26], where the basic equations are in Appendix A, the inelastic 
cross section is 



J d 2 b [1 - e" 2 ^] , (11) 
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FIG. 1: Minimum transverse momentum (p™ in ) required for the collision to qualify as a hard 
scattering, sibyll 1.7 had a constant minimum at v5 GeV, whereas this has been modified to 
change with yfs in version 2.1. 



where the eikonal is 

X(M) = Xhard(&, S) + Xsoft(M) = ~™hard(M) + -n soft (M) . (12) 

The number of soft interactions is denned analogously to the hard one n so f t (6, s) = 
c"soft(<s). The soft part of the eikonal is discussed in the following subsection. 
The hard part of the eikonal is interpreted as having a probability of exp[— n^dib, s)] for 
no minijet production at energy s and impact parameter b. Equation (11) can be reorganized 



as 
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(14) 



is the cross section for production of N pairs of minijet. This interpretation follows from 
the Abramovski-Gribov-Kancheli [27] cutting rules, where is the term with exactly N 
cut parton ladders, summed over all uncut ladders [28]. The probability distribution for 
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FIG. 2: Minijet production, with the energy fraction carried by the minijets (left panel) and the 
average number of minijets produced (right panel) over a range of yfs. 



obtaining N minijet pairs is 



(15) 



-*jV — ! 

^incl 

and the mean number of minijet pairs produced per interaction is 

oo 

(N) = Y. NP * = • ( 16 ) 

The contribution of minijets to the overall particle production for a p-p collision is shown 
in Fig. 2, where the energy fraction carried by the minijets and mean number of minijets 
produced are shown as a function of cm. energy. Minijets start becoming important at 
y/s « 1000 GeV. 

Each minijet pair is treated as two strings stretched between two gluons. In order to 
fragment the jets, a q-q pair is generated at each end, and a leading particle at each end is 
created. Then the string fragments in the standard way. The fraction of energy going into 
the minijet from each hadron 1, 2 {x\ an x 2 ) is obtained by selecting x from the effective 
parton density function [29] 

f(x) = g{x) + ^ [q(x) + q(x)} . (17) 

The code uses the approximation that for small transfer momentum the cross sections for 
g-g, q-g and q-q scattering are proportional to t~ 2 and are in ratio 1 : 4/9 : (4/9) 2 . The 
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transverse momentum follows 



da 1 
— oc — , 
dt t 2 



(18) 



where t is the four-momentum transferred squared Mandelstam variable, and i > (p™ in ) • 



We emphasize that the full parton structure functions of the u, d, s, c quarks and gluon are 
used for the calculation of the hard cross section. The above approximation is made only 
when sampling partonic final states and the parton density is parametrized in this simple 
way by adding quarks and gluons with the approximate weights. 

C. Soft interactions 

SIBYLL 1.7 has a very simple, energy-independent form of soft contribution to the eikonal 
Xsoft = \CA(b), having an impact parameter profile function identical to that of the hard 
counterpart, and C = 123 GeV -2 is chosen to reproduce the low energy inelastic cross section 
of 32 mb. Only one soft interaction is permitted, and the hard-soft interaction division was 
energy-independent at pt = V§ GeV. The energy left over after the production of minijets 
was shared by two strings connecting the valence quarks of the projectile and target. 

In version 2.1, the energy-dependent p™ m (s) allows a larger range of phase space for soft 
interactions. The eikonal form of Xsoft = ^A so{t (b)a so f t (s) is kept. We adopt some aspects of 
Regge theory in order to accommodate multiple soft interactions. Inspired by Regge theory 
[30] the energy dependence of cr so ft is taken as sum of two power laws, one term for Pomeron 
exchange and another term for Reggeon exchange [31] 



The index e for Reggeon exchange at low energy is expected to be very similar to the one 
found in fits by Donnachie and Landshoff [31]. The parameter A e g, in contrast, depends on 
the subdivision of the Pomeron term into soft and hard contributions and is hence a function 
of the transverse momentum cutoff (7). Here we implicitly assume that minijets form the 
hard part of the Pomeron [32]. 

The parameters X, Y and e and A e g are determined by fitting the measured total, elastic 
and inelastic cross sections for p-p and p-p interactions. Based on the GRV parton densities 
[20, 21] e « 0.4 and A off « 0.0245 are found. 




(19) 
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FIG. 3: Geometric interpretation of a hard interaction (left) and soft interaction (right). 



In the following we describe soft interaction by extrapolating the picture of hard interac- 
tions into the domain of low momentum transfer, in which this picture cannot be justified 
within perturbative QCD. Only comparison of the corresponding model predictions with 
data can later prove this description as useful. 

While hard interactions are approximately pointlike in character, soft interactions involve 
a larger transverse interaction area. The low px of the partons in a soft interaction spreads 
out the interaction area (from the uncertainty principle ApAb ~ 1), as opposed to a high 
event which localizes the collision to a small region. The geometry of hard and soft inter- 
actions in impact parameter space is schematically shown in Fig. 3. The energy- dependent 
increase of the interaction region of soft interactions has to be taken into account in the soft 
profile function of the hadron-hadron collision 

A y f (s, b) = j d 2 ^ d 2 b 2 d 2 h 3 A y (b t ) A z (b 2 ) A soit (s, b 3 ) 5 (2) (bx - b 2 + b 3 - b) , (20) 

where the profile function A y / Z is analogous to the hard interaction case. The fuzzy area of 
the right hand diagram in Fig. 3 is represented by A soft (s,b3), which is parametrized as a 
Gaussian with a energy-dependent width 



A^\sM) = , = n exp 



(21) 



with 



(22) 



B s (s) = B + a'(0)lnW . 

In the limit of B s (s) — > 0, A soft (s, bs) becomes a delta function and Eq. (20) becomes equal 
to the hard profile function. 

10 



In order to calculate Eq. (20), an exponential form factor which corresponds to a Gaussian 
shape in transverse space is used for the proton or meson profile function A y (b) as a first 
estimate, using data to fit the parameters. For a proton-proton collision, this yields 



A(2B P + B a {s)) 



(23) 



™ v ' ' Att(2B p + B 8 (s)) 
where B p characterizes the transverse size of a proton and is fitted to data. 

For generating the string configurations in inelastic events, the number of soft (N s ) and 
hard (Nh) interactions is sampled from (see also [33]) 

Ns ' Nh ~ J N,\ N h \ ' [ ' 

with the inelastic cross section given by 

Cincl = °~N s ,N h ■ (25) 

N s +N h >l 

The probability distribution Eq. (24) is tabulated during initialization of sibyll and later 
used to draw event configurations. 

In analogy to hard interactions, soft interactions are simulated by a pair of gluons which 
are fragmented the same way as a minijet pair. Thus, there is one valence quark string pair 
and n s — 1 gluon pairs. The only difference is the lower transverse momentum, see Eq. (4), 
and the distribution of the momentum fractions of these gluons, which are sampled from 
an 1/x distribution. This distribution corresponds to the one expected for a scenario of 
saturated gluon density in central collisions. A minimum mass of m so ft = 1 GeV is required 
for strings between sea quarks to regularize the singular part of the distribution and to 
ensure applicability of string fragmentation. One of the multiple soft interactions always 
involves the valence quarks and the momentum fractions are then sampled from Eq. (1). 

Implementing multiple soft interactions affects the model predictions at intermediate 
energy, which can be seen, for example, in the inelastic and total cross sections between 
yfe = 50 - 900 GeV in Fig. 4. 



D. Diffraction dissociation 



Diffraction is a collision where there are no quantum numbers exchanged between the 
colliding particles. A characteristic feature is a large rapidity gap in the final state. Unfortu- 
nately diffraction physics is not satisfactorily understood even on the level of phenomenology. 

11 



A comprehensive description of diffraction can be found in Ref. [34]. We give only a brief 
description of the diffraction model used in the new version and put the relevant equations 
in Appendix C. 

In version 1.7, diffraction was considered part of the inelastic, no-minijet event but was 
not otherwise included within the physics framework. The cross sections a^is were simply 9% 
each for forward and backward diffraction and 4% for double diffraction of the <Ti ne i (including 
minijet production) at 30 GeV, and was assumed to increase with energy as adis oc ln(s), 
with the diffractive event probability P^m = Cdifr/ow- With this treatment, however, cr diff 
becomes larger than the cross section for inelastic no-minijet events at high energy. Also, 
the minijet cross section is a fit to the total inelastic cross section, which includes diffractive 
events as well. This resulted in an underestimation of minijet production. 

SIBYLL 2.1 uses the two-channel eikonal model to incorporate diffraction into the eikonals 
for the low-mass diffraction dissociation. The procedure is similar to the Good- Walker 
model [14] with only a few assumptions made; see also Ref. [35]. Only two states are 
distinguished, a nonexcited state and a generic diffractively excited state (denoted with a 
*) which is a superposition of different low-mass diffractive final states. The cross sections 
and relevant equations are given in Appendix C. The high-mass diffraction dissociation is 
calculated by extracting the corresponding eikonal function from the soft and hard eikonal 
parts. Including the excited state of projectile and target in the eikonal formalism gave 
considerable improvement to the multiplicity distribution (see. Sec. Ill A). 

Diffraction dissociation is treated with a strict kinematic cutoff M^/s < 0.1, which follows 
from considerations on coherence and diffractive particle production [36]. The net effect is 
that the quasielastically scattered protons do not lose more than ~ 20 % at maximum 
in diffraction dissociation at this energy. The diffractively dissociating particle undergoes 
a phase-space decay if the mass of the excited system is very low. For higher masses, 
diffracted particles are divided into two valence components of quark-diquark or quark- 
antiquark that are connected by a color string which subsequently fragments. The string 
carries the diffracted particle's momentum and quantum numbers and does not create extra 
Pt- The one-string decay threshold is set to AM = 0.7 GeV, where AM is the mass 
difference of the incoming particle and the excited state of it. 

Diffractively excited states of a mass of more than 10 GeV are considered as being pro- 
duced by a Pomeron-hadron interaction. The decay of these states is described with multiple 
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soft and hard interactions by generating a tr-p interaction at = AM, as motivated by 
data from the UA4 Collaboration [37]. 



E. Nucleus interactions 

The physics framework for the hadron-nucleus and nucleus-nucleus interactions remains 
the same in the new version. The difference in the cross sections between the versions comes 
from the improvements made to the hadron-hadron interactions. Detailed descriptions can 
be found in Ref. [1] for the hadron-nucleus interaction and Ref. [12] for the nucleus-nucleus 
interaction. Here, we summarize the basic concepts for the sake of completeness. 

The interaction length is calculated from the production cross section. Elastic and 
quasielastic interactions, where no new particles are produced, do not contribute to the 
air shower development and are not considered. The production cross section, i.e. contribu- 
tion to particle creation, for a hadron-nucleus (hA) interaction where A is the mass number 
of the nucleus is 

„hA „hA „hA JxA I 

with of Q ^, a^ A , and being the total, elastic, and quasielastic cross sections, respectively. 
They are calculated within the Glauber model [11] from the p-p, n-p, and K-p cross sections. 
The inelastic and total cross sections of p-p, p-aii and 7r-air, ir-p collision are shown in Fig. 4. 
The minijet cross sections are also shown. 

In a hadron-nucleus interaction, the number of target nucleons directly participating in 
the interaction, also known as wounded nucleons, is determined from the production cross 
section. The mean number (N w ) of wounded nucleon per interaction is given by standard 
Glauber theory. In analogy with Eqs. (15) and (16), one finds 

i A hp 

"prod Nw "prod 

where gn w is the cross section for interaction with N w nucleons and Gh v is the hadron-nucleon 
cross section. 

The string model is applied to the fragmentation of the partonic system as well. The 
target nucleus is seen as N w pairs of valence q-qq, and the projectile hadron is viewed as one 
valence q-q or q-qq pair and N w — 1 sea q-q pairs. The N w color-connected partons undergo 
string fragmentation. Most of the energy is carried by the valence pair string, and the 
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FIG. 4: The total and inelastic cross section of p-air, p-p (left panel) and 7r-air, ir-p (right panel) 
collision for version 2.1 (red solid lines) and 1.7 (blue dotted lines). The minijet cross section is 
also shown in dashed lines, where the higher one is version 2.1. 

sea strings contribute to giving the proper multiplicity in the target region. This happens 
because there are more valence quarks on the nucleus side for the sea quarks to couple 
to. Intranuclear interactions and Fermi momentum inside the target nucleus are ignored 
in SIBYLL. Intranuclear interactions are greatly suppressed at high energy due to the time 
needed for a hadron to form as an independent object (formation time). At high energies, 
minijets are added to the collision, by generating them for each of the N w wounded nucleons 
in the same way as in a hadron-hadron interaction. 

The nucleus-nucleus interaction is treated with the semisuperposition model [12], which 
is between the simple superposition model and full Glauber theory [11]. The superposition 
model treats each nucleon of the projectile independently and as a consequence the interac- 
tion lengths of the nucleons have an exponential distribution based on the hadron-nucleus 
cross section. In reality, the nucleus interaction length is very small and a nucleus will inter- 
act quickly in the atmosphere. In the semisuperposition model, the number of interacting 
nucleons in the projectile for each nucleus interaction is determined from Glauber theory, 
where the remaining spectator nucleons fragment into lighter nuclei. Though the interac- 
tion and fragmentation is treated as a nucleon-nucleus interaction, the distribution of these 
nucleon-subshowers reflects correctly the nucleus-air cross section. 
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III. COMPARISON WITH EXPERIMENTAL DATA 



Both fixed target and collider experiments give valuable guidance in modeling hadronic 
interactions. Fixed target experiments provide data for the forward region which are most 
relevant to cosmic ray interactions, but the energies are relatively low, ~ several hundred 
GeV. Collider experiments can probe higher energies (£i ab ~ 10 6 GeV) but most of the 
information is collected in the central region. Some collider experiments such as HI and 
ZEUS are able to detect forward events [38]. The anticipated experiments LHCf [39] and 
TOTEM [40] in the LHC at CERN are expected to collect information in the forward region 
at an energy equivalent to cosmic rays of E\ ah ~ 10 8 GeV. 

A. Charged particle multiplicities 

Pions are the most numerous particles, followed by kaons and baryons. There is overall 
good agreement with experimental data in the forward region at low energies. The difference 
of charged particle production between the two versions is due to the improved treatment 
of multiple soft interactions, usage of GRV parton densities, and a consistent inclusion of 
diffraction dissociation, which also leads to more minijet production. These improvements 
give a better agreement with data for version 2.1, especially in the central region. 

The NA49 experiment measured the rapidity y and Feynman xf distribution of charged 
particles for p-p [41] and p-C [42] collisions at E^b = 158 GeV. Figure 5 shows the sibyll 
results compared to the data for ir + and n~, which became available only after the event 
generator had been released. Good agreement between model predictions and data is found. 
The excess of 7r + over tt~ is due to the flavor content of the proton (uud). In version 2.1, this 
discrepancy is stronger and more particles are produced in the central region which reflects 
the changes made to the soft interaction. However, the difference between the two versions 
is small in this respect. 

Fixed target experiments at FNAL used tt + , K + , p as projectiles and p, C for targets. 
The inclusive cross section Ed 3 a/dp 3 for each charged particle species has been measured at 
Eiab = 100 GeV at a given px [43]. The results for ti + ix~ production at pr = 0.3 GeV/c are 
plotted in Fig. 6. The inclusive cross section of version 2.1 is slightly lower than version 1.7, 
and gives an overall better agreement. For the pion projectile, 7r + s are overproduced while 
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FIG. 5: The Feynman x (xf) and rapidity (y) distribution of pions plotted against NA49 result of 
p-p [41] and p-C [42] collision at E\ a \, = 158 GeV. Version 2.1 (1.7) results are shown in red solid 
(blue dotted) lines. The left (right) panels show the production of ir + (vr~). The upper panels 
show the xp distribution, where the p-C collision results are multiplied by factor 10 in order to 
show both interactions on the same plot. The lower panels show the y distribution: the upper 
(lower) set of lines and data points are from the p-C (p-p) collision. 

7r~s are underproduced in the forward region, indicating the role of the valence components 
in determining the leading particle. The kaon projectile shows an overall good agreement 
for version 2.1. For the proton projectile, there is an overproduction of n~s compared to 
data, especially in the forward region. 

Pseudorapidity (rj) distributions of charged particles from collider experiments are com- 
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FIG. 6: Inclusive cross sections of n + -p,C (left panels), K + -p,C (center panels), p-p,C (right 
panels) collisions for tt + (upper panels) and ir~ (lower panels) production, where the events have 
Pt = 0.3 GeV. Collision energy is E^b = 100 GeV, and are compared with results from a FNAL 
fixed target experiment [43]. Version 2.1 (1.7) results are shown in red solid (blue dotted) lines. 

pared with sibyll in Fig. 7. It shows the t] distribution of charged particles from p-p 
collisions at E c . m . = 1800 GeV (CDF [44]), 630 GeV (P238 [45]), 200 GeV (UA5 [46]) and 
53 GeV (UA5 [47]). The improvements made to version 2.1 most prominently show in the 
central region. The role of the minijets and soft interactions is visible in the central region, 
where version 1.7 lacks secondary particles especially as the energy increases, while hav- 
ing more particles in the peripheral region. This trait can be seen at low energies in the 
pp — > 7r + , 7T~ figures in Fig. 6. Version 2.1 gives an excellent description of P238 data and 
tends to slightly overestimate the particles at low energies. It should be noted that the r\ 
range and trigger condition for 53 GeV is different than for higher energies at UA5. The two 
versions are similar for events with large \r)\ beyond the scope of current collider detector 
measurements. 

The distributions of charged particle multiplicity at UA5 [48] also give information at 
higher energies. Figure 8 shows the distribution of charged particle multiplicity for p-p 
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FIG. 7: Pseudorapidity (77) distribution of p-p collision, for various cm. energies, for version 1.7 
(left panel) and 2.1 (right panel) Data points from top are E c , m , = 1800 GeV [44], 630 GeV [45], 
200 GeV [46], 53 GeV [47]. Note the deficit of charged particles in version 1.7 at high energies in 
the central region. 




FIG. 8: Probability of charged particle multiplicity (n c h) distribution at UA5 detector with E c _ m , = 
900 GeV for three different rj ranges; \rj\ < 0.5 (xO.Ol), |r/| < 1.5, |r/| < 3 (xlOO) [48] for version 
2.1 (red solid lines) and 1.7 (blue dotted lines). The two sets of line are for having a stable and 
decaying K®. 
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collision at E c , m , = 900 GeV, at three different r] ranges \rj\ < 3.0, 1.5, 0.5, where the 
results for \r]\ < 3.0 have been multiplied by 100 and \r]\ < 0.5 by 0.01 for clarity while 
plotting. The particle K® has a very short lifetime and its decay produces charged particles 
that have a non-negligible effect, especially at high multiplities. The treatment of K® is 
considered as one of the uncertainties in the interpretation of the experimental data. Both 
stable and unstable cases are plotted, which can be considered as an error band. The 
improvements made in soft interaction and diffraction in version 2.1 are evident in the wider 
distribution of n c h as well as in the increase in multiplicity. An underestimation of the cross 
section for double diffraction dissociation in sibyll is probably the reason for the lack of 
low-multiplicity events satisfying the UA5 trigger condition. 

B. Leading particle 

The leading particle from the fragmentation carries a significant fraction of the total 
energy and becomes the primary particle in the next interaction of the air shower. The 
elasticity K = Ei cad /E pro j, the fraction of the leading particle with respect to the collision 
energy, of a collision affects the multiplicity as well as the speed of shower development in 
the atmosphere. Thus it is important to get a correct description of the behavior of the 
leading particles. 

The NAL bubble chamber experiment has data for p-p interactions at = 102, 205, 
303, 405 GeV and measured the xp of the leading proton [49]. Figure 9 shows the sibyll 
results plotted against the NAL data. The sharp dip at xp ~ 0.9 for the old version indicates 
the abrupt onset of diffraction, which is softened for version 2.1. It is not a smooth turn-on 
however, with the switch-on evident from the small step around xp = 0.8. 

The ZEUS detector measured the leading proton [50] with small transverse momentum. 
ZEUS collided positrons with protons, where the proton energy was 820 GeV during run 
I. The cm. collision energy is about 300 GeV. A virtual photon emitted from the positron 
interacts with the proton. The ZEUS Collaboration has confirmed that the initial projectile 
gives little effect [51] on the outcome. As sibyll cannot have a photon or positron projectile, 
we simulated a p-p collision at a slightly lower energy of -E c . m . = 210 GeV and used events 
from one hemisphere, i.e. events with p c z m ' > 0. Figure 10 shows the leading protons of the 
sibyll results plotted against the ZEUS data. They are plotted as a function of xi a b = E/E p , 
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FIG. 9: Distribution of leading protons as a function of xf compared to the p-p collision data from 
the NAL bubble chamber [49]. Version 2.1 (1.7) results are shown in red solid (blue dotted) lines. 
The results and datapoints have been multiplied by factor of tens (405 GeV by 1000, 303 GeV by 
100, 205 GeV by 10, 102 GeV by 1). 

the energy of the proton or neutron divided by the colliding proton energy in the lab frame, 
which is essentially the elasticity. The leading proton displays similar behavior to that of the 
NAL bubble chamber. Again, the better diffraction treatment is evident around xi a b = 0.9. 
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FIG. 10: Leading proton distribution of SIBYLL compared to the low-Q 2 [beam pipe calorimeter 
(BPC)] and higher-Q 2 [deep inelastic scattering (DIS)] ZEUS data [50, 51]. Version 2.1 (1.7) 
outcome is shown in red solid (blue dotted) lines. The transverse momentum of the protons are 
limited to p\ < 0.5 GeV 2 . The cm. collision energy of sibyll is 210 GeV. 

C. Strange particle production 

A FNAL fixed target experiment measured the production of very forward strange par- 
ticles produced in p-Be collisions at = 300 GeV [52]. The inclusive cross section of A 
and K® have been plotted for angles in the range 9 = 0.25 — 8.8 mrad in Fig. 11. As these 
are very forward direction measurements, particles with large p z are likely to be from par- 
ticle production associated with the projectile, and particles with small p z are from central 
production. Both sibyll versions give agreements in the forward direction, with a tendency 
to slightly overproduce high-p 2 particles and underproduce low-p^ particles. Strange particle 
production directly affects production of high energy muons and neutrinos. 

IV. AIR SHOWER PERFORMANCE 

The development of an air shower depends on a number of factors, some of which are 
the production cross section, inelasticity, and multiplicity. For the description of the air 
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FIG. 11: Inclusive cross sections of A (left panel) and K® (right panel) production from p-Be 
collision at = 300 GeV for various angles [52]. Version 2.1 (1.7) results are shown in red solid 
(blue dotted) lines. From the top, the sets of lines and data are for angles 0.25, 1.5, 2.9, 5.0, 6.9, 
8.8 mrad. The results and data have been multiplied by factors of tens to show them in one figure 
(0.25 mrad by 10 5 , 1.5 mrad by 10 4 , 2.9 mrad by 10 3 , 5.0 mrad by 100, 6.9 mrad by 10, 8.8 mrad 

by l). 

shower development the hadronic interaction model has to describe correctly the particle 
interactions in a wide range of energies. Observables such as depth of shower maximum X max , 
electron number N e and muon number at ground will depend on the characteristics of 
hadronic interactions 1 . We briefly summarize how air showers are affected by those three 
parameters. 

Increasing the cross section will cause the shower to start earlier in the atmosphere, re- 
sulting in a smaller X max as well as a smaller fluctuation. The number of electrons measured 
is highly dependent on the position of the shower maximum; the closer to X max the larger 

1 Extensive studies have been carried out in Rcf. [53], where cross sections, elasticity, and multiplicity are 
varied to see the effect on observables such as X max , N e and N^. 
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FIG. 12: The mean X max of protons and iron nuclei as a function of primary energy at zenith 
angle 6 = 45°. In the case of sibyll (red solid lines), QGSJET II-3 (orange dot-dashed lines), and 
epos 1.61 (black dashed lines), the data are taken from the hybrid simulator CONEX [54, 55] with 
statistical uncertainly of typically 1—3 gem -2 . For SIBYLL 1.6 (blue dotted lines), 500 events were 
simulated with Corsika 5.62 [56] (from Ref. [57]) with statistical uncertainty of 1 — 4 gem" 2 . 

the N e . As most muons are produced from decay of pions and kaons, is expected to 
remain stable. 

An increase in the mean multiplicity lessens the energy per particle of the secondaries, 
which results in a quicker development of the shower with smaller fluctuations. The in- 
creased multiplicity is expected to increase N^. The number of electrons are most numerous 
at A max and decreases away from it. Hence, a quicker development results in a larger dis- 
tance between A max and ground which will contribute to decreasing the number of electrons 
detected on ground. 

A large elasticity gives a larger fraction of the total energy to the leading particle, and 
slows the shower development as well as giving a smaller multiplicity. The A max will be 
closer to the ground and will be smaller. The closer proximity of A max to the ground 
causes a larger N e despite the smaller multiplicity [53]. 
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The new version has the shower developing more quickly than the old version. The two 
versions are compared with qgsjet II-3 [58] and EPOS 1.6 [59]. The shower maximum is 
shown in Fig. 12. sibyll 2.1 has an increased cross section and larger multiplicity which 
results in a smaller X max overall but maintains the same shape. QGSJET II-3 has a very 
large multiplicity at high energies, resulting in smaller shower maximum. 

The muon number is an important indicator for cosmic ray composition studies, as show- 
ers of heavier nuclei contain more muons than that of protons. Version 2.1 produces more 
muons than version 1.7, which cause the new version of sibyll to extract a lighter cosmic 
ray composition from experimental data than the old one. However, both versions produce 
fewer low energy muons than other models such as QGSJET II-3 or EPOS 1.6. The muon 
number at sea level is sensitive to zenith angle. When the atmosphere depth increases, 
muon production and muon energy loss and decay compete against each other. For low 
muon energies, the decay process is dominant and V M decreases as increasing zenith angle. 
At sufficiently high energies above ~ 10 GeV, most muons do not decay and increases. 
Figure 13 shows the energy dependence of the average number of muons normalized to the 
primary energy for the two versions of sibyll. The average number of muons at sea level 
((Nfj)) is plotted with energies above E^ r =0.3, 1, 3, 10 and 30 GeV in proton-initiated 
showers at zenith angle = 0°. 

V. SHORTCOMINGS AND FUTURE IMPROVEMENTS 

Though many features have been improved, sibyll 2.1 is by no means a complete model. 
We list the shortcomings of the current version and describe how we plan to improve the 
model. 

• The current nucleus-nucleus collision uses the semisuperposition model. Implementing 
the full Glauber model will give a more accurate description. 

• Antibaryon production is not satisfactorily described. There is not enough produced 
and the distribution of antiprotons in the central region is incorrect compared with 
data [61]. The overall normalization can be improved by increasing the diquark produc- 
tion fraction. The current method of fragmentation suppresses antibaryon formation 
together with other particles in the non-end of strings. 
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FIG. 13: Average number of muons at sea level (Nu), obtained in proton showers with zenith 
angle 9 = 0°. Each energy represents 5000 showers simulated with the hybrid method. The solid 
(dotted) line represents the values obtained with sibyll 2.1 (sibyll 1.7). Panels (a)-(c) show the 
average number of muons with energy above 30, 3 and 0.3 GeV respectively. Figure is adapted 
from Fig. 14 of Ref. [60]. 

• The currently used energy-dependent transverse momentum cutoff is independent of 
the relevant gluon density of the interaction with the target nucleus, which varies with 
the impact parameter of the collision. A new energy- and impact parameter-dependent 
cutoff to pt to prevent parton density saturation would improve the modeling. Refer- 
ence [62] is an example of an attempt to constrain minijet formation. In addition the 
profile functions could also be more refined. 

• A consistent treatment of coherent and incoherent diffraction dissociation in hadron- 
nucleus and nucleus-nucleus interactions is required. This can be achieved by using 
a two-channel model in the Glauber calculation similar to the one presented here for 
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p-p interactions. 

• Include charm quark. The current model has u, d s quarks and gluons. This im- 
provement will be relevant more to muon and neutrino detectors than cosmic ray 
observations. A version containing charm will be released soon. 

VI. CONCLUSION 

In this paper, we described the overall model of sibyll together with the changes made 
in the version 2.1, and listed the shortcomings and possible ways of improving the model. 
The 2.1 version still keeps the DPM-minijet structure but with modifications and additions. 
Results from HERA suggest a steeper parton density for gluons at low x, and the GRV 
parton densities replaced the EHLQ parton densities which resulted in an increased QCD 
cross section. Concepts from Regge theory are used to allow multiple soft interactions. The 
energy-dependent transverse momentum cutoff for ensuring perturbative QCD to be valid 
is a better discriminator between soft and hard interactions than the previously applied 
energy-independent cutoff. The two-channel eikonal model is used to describe diffraction, 
with better if not complete success. The physics framework for the hadron-nucleus and 
nucleus-nucleus interactions remains unchanged. 

These improvements produce more particles with a wider distribution in momentum 
space, as are evident when comparing with experimental data. Both versions give a good 
fit to the rapidity and Feynman x distribution for fixed target experiments such as NA49. 
However the changes made to the new version are very evident in the central region, in the 
pseudorapidity and the overall charged particle multiplicity distribution. As a consequence, 
air showers described with sibyll 2.1 develop quicker, with smaller shower maximum and 
larger muon number than version 1.7. 

SIBYLL 2.1 is by no means the final answer. There are still shortcomings, such as un- 
satisfactory description of antibaryon production, and further room for improvements in 
preventing parton density saturation by giving an energy-dependent cutoff to the transverse 
momentum, giving a better diffraction treatment and have a full Glauber model for nucleus- 
nucleus interaction. A version including charm quark production is due to be released 
shortly. 

With all its advantages and shortcomings, the interaction model sibyll 2.1 is able to 
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successfully reproduce and predict many experimental data. We believe there that use of 
SIBYLL in comparison with other models would be valuable in analysis of cosmic ray air 
shower data. 
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APPENDIX A: AMPLITUDE CONVENTIONS 

The conventions and parameters are adopted from Ref. [26]. The Mandelstam variables 
s, t, u are used, and k is the cm. momentum. The cm. scattering angle 9 is related to t by 

t = -Ak 2 sm 2 (9/2) , (Al) 

and a new parameter is defined — q 2 = t. The scattering amplitude in the cm. frame 
[fc.m.( s i t)] an d the Lorentz invariant scattering amplitude [^(s,t)] are related by 

Jl{s, t) = -8tt / c . m . (s, t) . (A2) 

The scattering amplitude f c .m.( s ^) an d impact parameter function a(s,h) are Fourier 
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transforms of each other via 



/ c .m.(M) = ^J d 2 be l * h a(s,b) , 

a(s,b) = -^J d 2 qe^ h f c . m .(s,t) . 



(A3) 
(A4) 



With this scattering amplitude, the differential elastic cross section can be expressed as 
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where the last relation used the optical theorem. The elastic slope parameter B e \ is defined 
from an approximation of the elastic scattering cross section in small t region as 
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Using Eqs. (A5a) and (A5b), 
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where p is the real to imaginary ratio of / c . m .(0). 

When using eikonals, the cross sections follow the usual convention: 

fftot = 2vr J db 2 (l - e"*^) 



0"el 
Cinel : 
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db 2 (l - e~ 2 ^) . 
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We neglect the real part of the elastic scattering amplitude for calculating the eikonal func- 
tions. 



APPENDIX B: HARD INTERACTION MINUET CROSS SECTION AND PRO- 
FILE FUNCTIONS. 

Each proton or meson is characterized by a (transverse) density profile function A z (b), 
where z can be p, ir, K. The probability for the two partons in particles y and z to collide 
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is found by integrating over all possible impact parameters k>i and b2 for a given impact 
parameter b of the collision 



b) = J d 2 b a d 2 b 2 A y (v y , bi) A z (u g , b 2 ) 5< 2 )(b a - b 2 - b) . (Bl) 
The profile function of a proton is given by 



(2vr) 2 7 V ^/ 

= i/„|i/ p b|lfi(|i/pb|) (B2) 

with i/ 2 0.71 (GeV/c) 2 . The profile function of a meson is 

Am{Vm: b) . _i_ | dV + ^ , (B3) 

where z/ m and 77 are adjustable parameters. For pions, v 2 « 0.54 (GeV/c) 2 and 77 ~ 0. Hence 
the profile function for a pp interaction is 



4r%> b ) = /rfv^db-b'D^db'D 



^ 1 



12tt8 

and for a /J7r interaction 



J- - (u p b) 3 K 3 (u p b) , (B4) 
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pn v p ' " ' 2tt (1 - C) 
where £ = (vp/v^) 2 - The profile functions are normalized to 



(B5) 



J d 2 bA^ d (u p ,u n ,b) = J d 2 bA y (v m ,b) = 1 (B6) 
APPENDIX C: DIFFRACTION DISSOCIATION 

In low-mass diffraction, the two-channel eikonal model is used [15, 35]. Only two states 
are considered here; a nondiffractive state and a generic diffract ive state denoted *. The 
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diffractive scattering of particles Y and Z can be expressed in the following matrix elements 



(YZ\^ int \YZ) 
(YZ\^' mt \Y*Z) 
{Y Z\JZ' ult \Y Z*) 
(Y 'Z\JT A \Y* 'Z*) 
(Y*Z\^' mt \Y*Z) 
(YZ*\^ int \YZ*) 
(Y*Z*\^ int \Y*Z*) 



^# Born 

fa ^ Born 

(3 z ^ Bom 

p Y p z ^ Boin 

(1 - 2a Y ) J( Bom 

(1 - 2a z )^ Boin 

(1 - 2a Y ) (1 - 2a z ) 



✓Born 



(Cla) 
(Clb) 
(Clc) 
(Cld) 
(Cle) 
(Clf) 
(Clg) 



The coefficients a and (3 may depend on energy. A matrix x(s, b) for the eikonal \ is 
introduced to calculate ' . The eikonal matrix is diagonalized, and the cross sections 
calculated. The hadronic states Y and Z are defined as 



\y,z) 






W 

The eikonal matrix reads 
/ 1 



Y*,Z) 



M 
i 
o 

W 



\Y,Z*) 





1 

w 



\Y*,Z*) 






VV 



(C2) 



\ 



f3z (1 - 2a Y ) 
(3y (1 - 2a z ) 
- 2a Y ) (1 - 2a z )J 



X (s,b) . (C3) 



/3y /5z 
1 - 2a y /3 y /3 z 

Pz (3y(3z 1 - 2a z 

Pz(l - 2a y ) /3 y (1 - 2a z ) 

After diagonalizing x(s,b), the cross sections can be calculated. The total cross section is 
given by 

<jf z = 2 J d 2 b (YZ\ (1 - e"* (s ' b) ) \YZ) 

= 2 1 d 2 h j2f^ n r z 1 ' n (-i) n - Ms,h)r 

J n = l 



IV. 



(C4) 



where 



J j 



el.n 



1 - ^)(1 - «i - 7,) n + (1 + ^)(1 - a,- + 7,)" 



7i 



7j 



(C5) 
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and 



3 = Y,Z . 



7, = V^M 2 
Consequently the elastic cross section reads 

„el / j2u \ y re\,nre\,n/ i \n— 1 

a YZ = db 2^Jy Jz (-1J -j 

n=l 

The cross section for single diffraction dissociation of particle Y follows from 

J d 2 b \(Y*Z\ (1 - e"* (b ' s) ) \YZ)\ 2 

[d 2 h f2f^' n fz ,n (-ir' Mh \ s))n 

J 1 n - 



(C6) 



(C7) 



SD,y 



n=l 



(C8) 



using 



/ 



diff,n 



2 7j 



(C9) 



Finally the expression for double diffraction dissociation is given by 

-x(b, s )) \YZ)\ 2 



a YZ 



J d 2 b\(Y*Z*\ (1 - 



d 2 b 



E/y iff ' n /z iff,n (-i; 



-i(x(M) n 



n=l 



III 



(CIO) 



Note that after carrying out the sum over n in Eqs. (C4), (C7), (C8), (CIO) the cross sections 
can be written as impact parameter integrals over a sum of exponentials. 

The parameter range for otj and (3j is limited by the unitarity constraint that all eikonal 
functions have to be non-negative 



i_ a ._ 7 .> 



ctj < 1/2 



f3j > 



(Cll) 



A good description of the data is found for a = 0.2 and f3 = 0.5. 

The partial cross sections for N s soft and hard interactions follow from 



J k=l 



[2A fcX soft(s, b)]^ [2X kXh ^(s, b)] N » 



exp {-2A fc (xsoft(s, b) + Xhard(s, b))} 

(C12) 
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with 



Ai 



_ Oty \ I a z 
1Y J V 1Z 

lY J \ 1Z 



A 



7y/ V 7z 



and 

Ai = (1 - a y - 7y)(l - a z - 7z) 

A 2 = (1 - a Y - 7*0(1 _ a z + 7z) 

A 3 = (1 - a Y + 7y)(! _ a z - Jz) 

A 4 = (1 - a Y + tv) (1 - «z + 7z) • ( C14 ) 

For high-mass diffraction, it is assumed that a constant fraction of each cut (soft or hard 
interaction) corresponds to an rapidity-gap final state. The corresponding cross section is 
written, see Eq. (C12) 



a 



SD 
hm 



5 (cr 1;0 + cr ,i) , (C15a) 



DD c-2 / i i \ i n2 SD,Z . a 2 SD,K /.~m r-i \ 

^hm = <5 (o-i,o + o-o,i + 0i,i) + Pya^' + (3 z a lm . (C15b) 

The factor 6 is estimated by comparing with HERA data: 10% of all deep inelastic scattering 
events at low x correspond to diffraction (5 ~ 0.1). The last two terms in Eq. (C15b) 
represent the cross section for low-mass - high-mass double diffraction dissociation. 
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